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We  present  a  stable  method  to  recursively  linearize  the  acoustic  inverse  scattering 
problem.  It  turns  out  that  the  iU-posedness  of  the  problem  can  be  beneficially  used 
to  solve  it.  It  means  that,  due  to  iU-posedness,  not  all  equations  in  the  nonlinear 
system  are  strongly  nonlinear,  and  that  when  solved  recursively  in  a  proper  order, 
they  can  be  reduced  to  a  collection  of  linear  problems. 

Our  method  requires  solution  of  a  series  of  forward  scattering  problems  with  as¬ 
cending  wave  numbers  (or  frequencies).  At  each  frequency,  a  linear  least-squares 
problem  is  solved  to  obtain  an  approximate  forward  model  which  produces  the 
prescribed  scattering  data.  The  robustness  of  the  procedure  is  demonstrated  by 
several  numerical  examples  in  the  inversion  of  the  Helmholtz  equation  in  two  di¬ 
mensions. 
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1  Introduction 


The  solution  of  the  inverse  scattering  problem  requires,  in  essence,  an  inversion 
of  a  nonlinear  mapping.  There  are  two  major  difficulties  associated  with  this 
nonlinear  problem  in  two  and  three  dimensions  :  ill-posedness  and  local  minima, 
neither  of  which  has  been  addressed  satisfactorily.  The  fact  alone  that  the  nonlin¬ 
ear  system  at  high  wave  number  is  extremely  oscillatory  and  therefore  possesses 
numerous  local  minima  renders  the  most  popular  approach — nonlinear  optimiza¬ 
tion  and  its  various  modifications — fundamentally  unreliable.  In  another  notable 
approach  to  attacking  the  nonlinear  problem,  the  original  boundary  value  prob¬ 
lem  is  reformulated  as  an  initial  value  problem.  The  resulting  method  (known 
under  several  names  such  as  layer-stripping,  trace  method)  avoids  the  problem 
of  local  minima,  only  to  be  plagued  by  ill-posedness.  There  are  still  other  ap¬ 
proaches  which  compromise  in  some  way  but  have  not  yet  been  able  to  squeeze 
through  the  flails  of  these  two  perils. 

We  present  a  stable  method  that  solves  the  inverse  scattering  problem.  It 
turns  out  that  the  ill-posedness  of  the  inverse  problem  can  be  beneficially  used 
to  solve  it.  It  means  that  not  all  equations  in  the  nonlinear  system  are  strongly 
nonlinear,  and  that  when  solved  recursively  in  a  proper  order,  they  can  be  reduced 
to  a  collection  of  linear  problems. 

Our  method  requires  solutions  of  a  series  of  forward  scattering  problems  with 
ascending  wave  numbers  (or  frequencies).  At  each  frequency,  a  linear  least- 
squares  problem  is  solved  to  obtain  an  approximate  index  of  refraction  which 
not  only  minimizes,  but  also  eliminates  the  discrepancy  between  the  prescribed 
and  the  predicted  scattering  data.  The  robustness  of  the  procedure  is  illustrated 
by  several  numerical  examples  in  the  inversion  of  the  Helmholtz  equation  in  two 
dimensions. 

The  plan  of  the  paper  is  as  follows:  in  Section  2  we  summerize  the  relevant 
analytical  apparatus,  in  Section  3  we  reformulate  the  ill-posedness  and  the  in¬ 
verse  scattering  problem,  and  present  the  inversion  method,  and  in  Section  4  we 
describe  the  numerical  implementation  of  the  algorithm.  The  robustness  of  the 
procedure  is  demonstrated  with  numerical  examples  in  Section  5. 


2  The  Mathematical  Preliminaries 

In  this  section  we  summarize  several  analytical  results  regarding  the  scattering 
problem,  the  solution  of  a  nonlinear  problem,  and  the  solution  of  linear  ordinary 
differential  equations  of  matrices. 
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2.1  The  Scattering  Problem 

The  subject  of  this  paper  is  the  inverse  scattering  problem  for  Helmholtz  equation 
in  two  dimensions 

^<j>{x,y)  +  k'^{l  +  q{x,  y))(i>{x,  y)  =  0.  (1) 

In  (1),  ^  is  a  real  number,  q  a  smooth  function,  with  q{x)  >  —1  for  all  x  ^  B?. 
We  will  be  referring  to  k  as  the  wave  number  or  frequency,  to  the  function  q  as 
the  scatterer  or  the  forward  model.  We  assume  that  the  support  of  g  is  a  disk 
D{w)  defined  by  the  formula 

D{w)  =  {{x,y)\x'^ +  y'^  ]  (2) 

for  some  n?  >  0.  We  will  be  considering  solutions  of  equation  (1)  of  the  form 

(i){x,y)  =  <f)o{x,y)  +  'il;{x,y),  (3) 

with  ^0  :  D{zo)  — >  C  a  solution  of  the  equation 

/^4>q{x,  y)  +  k‘^4)o{x,  y)  =  0  (4) 

(to  which  we  will  be  referring  as  the  homogeneous  Helmholtz  equation),  and 
tp  ■  B?  C  the  solution  of  the  equation 

A^(x,  y)  +  A:2(l  +  q{x,  y))^(x,  y)  =  -k‘^q{x,y)(pQ{x,  y),  (5) 

subject  to  the  outgoing  (Sommerfeld)  radiation  condition 

—  ikip^  —  0.  (6) 

We  will  be  referring  to  p  as  the  total  field,  to  <f>o  as  the  incident  field,  and  to  ip  as 
the  scattered  field.  Furthermore,  given  an  incident  field  po  we  will  be  referring  to 
the  determination  of  the  corresponding  scattered  field  as  the  (forward)  scattering 
problem. 

Remark  2.1  We  measure  the  size  of  the  scatterer  by  the  number  of  wavelengths 
through  the  longest  ray  tube.  More  precisely,  suppose  that  the  curve  £  C  D{w)  is 
the  longest  ray  tube,  and  ds  is  the  arc  length  element  on  1;  then  the  number  of 
wavelengths  across  I  is  given  by  the  formula  is 

^  ^  +  q{x,y)ds.  (7) 

When  the  medium  is  quite  homogeneous,  namely,  when  q  is  small,  and  thus  the 
ray  path  is  roughly  straight,  (7)  is  reduced  to 


2.2  The  Scattering  Matrix  and  Scattering  Data 

In  this  subsection,  we  reformulate  the  scattering  problem  introduced  in  the  pre¬ 
ceding  subsection  via  the  scattering  matrix.  For  a  more  complete  description  and 
analysis  of  the  scattering  matrix,  see  [11].  Denote  by  Jm  the  Bessel  function  of 
the  first  kind  of  order  m,  by  Ym  the  Neumann  function  (or  the  Bessel  function 
of  the  second  kind)  of  order  m,  and  by  Hm  Hankel  function  of  the  first  kind  of 
order  m,  so  that 

Hm  =  Jm+i-Ym-  (9) 

Given  a  positive  number  z,  we  will  denote  by  Xz  the  linear  space  of  all  two-sided 
complex  sequences  {am}  such  that  for  some  c  >  0 

\am'Jm{_^')\  C,  (19) 

for  all  integer  m.  We  will  denote  by  Yz  the  linear  space  of  all  two-sided  complex 
sequences  {^m}  such  that  for  some  c  >  0 


for  all  integer  m. 


(11) 


Remark  2.2  It  is  well-known  (see,  for  example,  [4],  [5],  [3])  that  once  the  pos¬ 
itive  number  v  is  greater  than  z,  the  Bessel  function  J„{z)  strictly  decays  as  v 
grows,  and  that  the  Neumann  function  Yu{z)  strictly  grows  as  v  grows.  Moreover, 
Ju{z)  becomes  very  small,  and  Yi,{z)  becomes  very  large,  when  v  reaches 


No{z)  =  z-l-0  (z^)  . 


Finally,  for  v  >  No{z), 


1 

\/2jrm 


(12) 


(13) 

(14) 


Remark  2.3  Formulae  (10)  and  (11)  can  be  rewritten  as 
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It  is  well-known  (see,  for  example,  [1],  [11]  for  details)  that  inside  disk  D{w),  the 
solution  <j)Q  of  the  homogeneous  Helmholtz  equation  (4)  can  be  expressed  as 

OO 

<^o(r,^)=  E  (Xr,-Jm{krye^\  (17) 

m—  —  oo 

with  the  sequence  a  =  {am}  €  Xkw,  and  that  outside  D{'czi),  the  scattered  field 
■0  can  be  expressed  as 

CO 

v»(r,^)=  ^rn-Hm{krye”^^.  (18) 

m=— OO 

with  the  sequence  /?  =  {^m]  €  Fto-  The  well-posedness  of  the  scattering  problem 
implies  a  linear  mapping  5^,^  :  Xk-w  ^  Ykw  defined  by  the  formula 

=  S^,k-a,  (19) 

for  every  a  =  {am}  G  Xkrz,  of  the  incident  field  (17)  and  ^  —  {^m}  €  Ykw  of  the 
corresponding  scattered  field  (18).  The  linear  mapping  is  refered  to  as  the 
scattering  matrix  corresponding  to  the  scatterer  q  in  D[w).  For  a  fixed  k  >  0,  the 
scattering  matrix  is  evidently  all  we  can  acquire  from  scattering  measurements 
outside  the  disk  D{zc).  The  set  of  matrices 

{  5^7, fc  1  0  <  fc  <  OO  }  (20) 

is  all  the  information  we  can  collect  from  real-frequency  measurements,  and  is 
defined  as  our  scattering  data. 

Definition  2.4  An  aperture  (of  acoustic  measurement)  is  an  area  in  the  square 
[0,2x1  X  [0, 2x];  the  full  aperture  is  the  entire  square.  A  Fourier  aperture  is  an 
area  in  the  (two-dimensional)  Fourier  space,  and  for  such  an  area  E,  a  function 
g  €  ii^(H)  is  said  to  be  specified  on  a  Fourier  aperture  E  if  its  Fourier  transform 
gm,n  is  known  for  all  points  {m,n)  G  E. 

Remark  2.5  The  knowledge  of  the  scattering  matrix  SrD,k  is  equivalent  to  that 
of  the  full-aperture  measurements  taken  outside  D{w),  that  is,  the  acquisition  of 
each  scattered  field  outside  D{w)  corresponding  to  every  possible  incident  field 
4>o;  see  (56),  (57)  for  further  details. 

Remark  2.6  Obviously,  for  a  given  positive  number  r,  the  above  analysis  and 
definition  can  be  applied  to  a  scatterer  located  on  the  disk  D(r).  In  particular, 
there  is  a  scattering  matrix,  denoted  by  ST,k)  corresponding  to  this  disk  scatterer. 

Remark  2.7  For  two  positive  numbers  ri  <  r2,  a  scatterer  in  the  disk  D{r-[)  can 
always  be  regarded  as  a  scatterer  in  the  larger  disk  D{r2).  Therefore,  there  are 
two  scattering  matrices,  Sri,k  o-'^^d  Sr2,k,  corresponding  to  the  same  scatterer  on 
the  two  concentric  disks.  It  is  easy  to  verify  that  the  two  matrices  are  identical. 
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2.3  The  Riccati  Matrix  Equation 


We  show  in  this  subsection  that  the  scattering  matrix  5c, fc  is  a  member  of  a  family 
of  scattering  matrices  governed  by  a  Riccati  equation.  For  a  more  complete 
discussion  of  the  Riccati  matrix  equation,  see  [11].  For  r  >  0,  following  the 
standard  procedure  of  invariant  embedding,  we  define  the  chopped  scatterer  qD{r) 
by  the  formula 


qD(r){p,0) 


q{p,e)  if/)<r, 
0  if  p  >  r. 


(21) 


We  will  denote  by  ST,k  •  ^kr  ^  ykr  the  scattering  matrix  (see  Remark  2.6) 
corresponding  to  the  chopped  scatterer.  In  other  words,  for  every  o:  €  Xkr  there 
exists  I3  €  Ykr  such  that 

^  =  Sr,k-<^  (22) 


where 

OO 

W)=  T,  (23) 

m=— OO 


is  the  scattered  field  outside  D{r)  corresponding  to  the  incident  field 


<i>oM  =  E  c^m-Mkpye^^^  (24) 

m=— OO 


upon  the  scatterer  qD(r)-  For  any  r  >  w, 

Sr,k  =  5c, fe,  (25) 

due  to  Remark  2.7.  At  r  =  0,  the  chopped  scatterer  qn^r)  is  identically  zero;  any 
incident  field  generates  no  scattered  field.  Thus, 

So,k  =  0.  (26) 

Now  we  need  the  following  notation  in  order  to  present  a  ordinary  differential 
equation  for  Sr,k-  Denote  by  J^,  ifz  the  diagonal  matrices 

=  diag{...,J-i{z),Jo{z),Ji{z),...},  (27) 

H;,  =  diag{...,H-i{z),Ho{z),Hi{z),...},  (28) 

and  by  F  the  Fourier  transform  converting  a  function  in  X^[0, 2x]  to  its  Fourier 
coefficients  in  Let 

Ur)  =  ^  9{r,«y’"‘d9,  (29) 

for  a  smooth  function  g  :  B?  C,  and  let  gr  be  the  matrix  whose  (m,  n)-th  entry 
is  defined  by  the  formula 

(Pr)m,n  —  dm—nij'y  (^^i) 
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Remark  2.8  Denoting  by  :  i^[0,2jr]  L^[0,2n]  the  diagonal  linear  magging 


(9r-fm  =  9(r,»)-m  (31) 

for  all  f  G  X^[0,27r]^  we  can  easily  verify  that 

Qr  =  F^gr^F~^.  (32) 

Obviously,  for  the  scatterer  function  qr  is  a  Toeplitz  matrix  whose  first  quadrant 
(entries  with  both  row  and  column  indices  positive)  is  of  the  form 


9o(r) 

9-1 (^) 

9-2(7') 

9-3(7’) 

9-4(7’) 

9-5(7-)  •  •  •  ■ 

9i(^) 

9o(r) 

9-1(7’) 

9-2(7’) 

9-3(7-) 

9-4(7-)  ••• 

92(r’) 

9i(^) 

90(7') 

9-1(7’) 

9-2(7’) 

9-3(7-)  ••• 

93(r) 

92(1-) 

91(7’) 

90(7’) 

9-1(7-) 

9-2(7-)  ••• 

94(r) 

93(1-) 

92(7-) 

91(7’) 

90(7-) 

9-1(7-) 

95(r) 

94(1’) 

93(r) 

92(7-) 

91(7-) 

90(7-)  ••• 

As  a  function  of  r,  the  scattering  matrix  Sr^k  satisfies  a  Riccati  equation,  see  [11] 
for  details. 

Lemma  2.9  For  any  A:  >  0  and  all  r  >  0^  the  scattering  matrix  Sr^k  •  Xkr  ^  Ykr 
is  a  solution  of  the  Riccati  equation 

dSr  k  o  ^ 

"  Sr,k^HkT)'qr\Hkr'Sr^k  +  Jkr)^  (34) 

In  the  reminder  of  this  paper,  we  assume  that  the  initial  value  problem  defined 
by  (34)  and  (26)  is  well-posed,  and  that  the  entries  of  the  matrix  Sr^k  depend  on 
(r,  k)  smoothly. 

Remark  2.10  Given  k  and  q^  the  initial  value  problem  of  the  Riccati  equation 
(see  (34)  CLud  (26))  will  be  solved  in  the  forward  modeling  of  the  inversion  algo¬ 
rithm  described  in  Section  3.6;  see  (146)^  (14'^)  f^'^  details. 

Remark  2.11  According  to  formulae  (10)  and  (11),  the  sequence  a  =  {ctm}  i'^ 
(22)  may  grow  as  fast  as  J^^{kr),  whereas  the  sequence  (3  =  {(3m]  in  (22)  must 
decay  at  least  as  fast  as  i7“^(A:r).  Therefore,  Remark  2.2  and  the  fact 

Pm  —  (5*r,A:)m,n’^7T,  (35) 
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implies  that 


=  0{HZ'{kr)),  (36) 

=  O(Mkr)),  (37) 

for  an  arbitrary  integer  n  and  large  integer  m  >  No{kr).  Thus,  the  entries 
of  Sr,k  whose  row  or  column  indices  are  greater  in  absolute  value  than  No{kr) 
are  essentially  zero  (see  Figure  1);  the  dimension  of  scattering  matrix  Sr^k  is 
effectually  finite:  it  is  essentially  a  square  matrix  of  dimension  2-No{kr). 


The  structure  of  the  first  quadrant  (entries  with  both  row  and  column  indices 
positive)  of  the  matrix  Sr,k  is  shown  in  Figure  1.  In  the  figure,  a  part  of  a  matrix 
is  labeled  zero  to  indicate  that  the  entries  in  that  part  of  the  matrix  are  essentially 
zero. 


Figure  1:  Scatterer  on  D{r)  and  Structure  of  Its  Scattering  Matrix  Sr,k- 

Remark  2.12  Given  an  operator  T  which  maps  from  one  linear  space  X  to 
another  Y ,  we  will  casually  use  the  term  virtually-null  space  of  T:  an  element 
f  E  X  belongs  to  the  virtually-null  space  if  is  small.  We  will  also  loosely 
use  the  phrase  virtual  rank  of  an  operator  in  accordance  with  the  term  virtually-null 
space;  the  virtual  rank  of  the  scattering  matrix  Sr,k,  for  instance,  is  no  greater 
than  2-No(r,k)  (see  Remark  2.11).  We  may  further  think  of  virtual  rank  as  a 
real  number. 

Remark  2.13  Given  positive  numbers  k  and  r,  according  to  Remark  2.2,  the 
diagonal  matrix  Jkr  is  essentially  a  square  matrix  of  dimension  2-No(kr).  There¬ 
fore,  up  to  0[(A:r)^],  the  matrix  Jkr  can  be  regarded  as  to  have  the  same  structure 

as  that  of  Sr, k  (see  Figure  1). 

\ 

Remark  2.14  Given  positive  numbers  k  and  r,  combining  (35)  with  formulae 
(10),  (11),  we  obtain 

{Sr,k-Hkr)n,m  =  O  [Hjfikr])  , 

(,F[kr'Sr,k')m,n  ~  i,Jm(,kr'j'j  , 
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(38) 

(39) 


for  an  arbitrary  integer  m  and  large  integers  n  >  No{kr).  In  other  words,  the 
matrix  SrxHkr  only  has  2*iVo(^^)  effectually  nonzero  rows — the  rest  of  the  rows 
whose  indices  are  greater  in  absolute  value  than  NQ{kr)  are  essentially  zero.  Sim¬ 
ilarly,  the  number  of  effectually  nonzero  columns  of  Hkr‘Sr,k  olso  2-A^o(^?")* 


No{kr) 


S-H 


Figure  2:  Structures  of  the  Matrices  Sr^k'Hkr-,  Hkr^Sr^k  * 

Figure  2  shows  the  structures  of  the  two  matrices  Sr,k'Hkr  and  Hkr'Sr^k-  Only  the 
first  quadrant  of  the  matrices  (entries  with  both  row  and  column  indices  positive) 
is  depicted  here.  In  the  figure,  a  section  of  a  matrix  is  labeled  zero  to  indicate 
that  the  entries  in  that  part  of  the  matrix  are  essentially  zero. 

Remark  2.15  Given  a  frequency  k  >  0,  the  forward  scattering  problem  (see 
Section  2,1  and  Remark  2.10)  defines  a  nonlinear  mapping  from  the  scatterer  q 
on  D{vj)  to  the  scattering  data  S^j^k-  The  initial  value  problem  defined  by  (Sf) 
and  (26)  turns  out  to  be  a  convenient  tool  for  analyzing  and  numerically  treating 
this  nonlinear  system.  According  to  Remark  2,11,  S^^k  only  has  (2'N^{kw)Y 
entries  essentially  nonzero;  therefore,  there  will  be  the  same  number  of  equations 
in  the  nonlinear  system. 

Finally  in  this  subsection,  we  show  some  basic  algebraic  properties  of  the 
matrix  Riccati  equation  and  the  scattering  matrix. 

Definition  2.16  Suppose  that  N  is  a  positive  integer  or  +oo.  For  a  square 
matrix  A  =  (  —N  <  i,j  <  N  )  of  dimension  2N  —  1,  define 

B  =  ^N<zJ<N),  (40) 

hi  =  (41) 

Clearly,  the  operator  (•  •  -y  defines  a  linear  mapping  from  (^(2-^v-i)x(2N-i)  itself. 
The  following  lemma  is  a  direct  consequence  of  Definition  2.16. 
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Lemma  2.17  Suppose  that  N  is  a  positive  integer  or  +oo.  Suppose  further  that 
A,B  are  two  matrices  and  I  is  the  identity  matrix  in  C^'^N-i)x(2N-i) _  finally, 
suppose  that  T  €  is  the  diagonal  matrix  T,-,i  =  (—!)*•  Then 


(A  •  BY 

II 

to 

A 

(42) 

(AT 

=  A, 

(43) 

r 

=  Y 

(44) 

=  T. 

(45) 

The  following  lemma  is  a  direct  consequence  of  Lemmas  2.17  and  2.9. 

Lemma  2.18  Suppose  that  r,  k  are  positive  numbers,  and  that  ST,k  is  the  scat¬ 
tering  matrix  corresponding  to  the  chopped  scatterer  qD(r)-  Suppose  further  that 
T  is  diagonal  matrix  of  infinite  dimension  defined  by  Tt,,-  =  (—1)*.  Then 

(JkrY  =  r-Jkr  =  Jkr-r, 

{HkrY  =  T-Hkr  =  Hkr-T, 

{q^Y  =  T-9.-T, 

{Sr,kY  =  Sr,k- 

2.4  The  Far  Field 

In  this  subsection,  we  first  define  the  far  field  of  a  scattered  field,  and  link  it  to  the 
scattering  matrix.  We  then  introduce  the  translation  operator  of  the  scattering 
matrix.  We  will  require  the  column- vector  Fg  and  the  row-vector  defined  by 


Fg  =  {e-*”^^m  =  ...,-2,-l,0,l,2,...  (50) 

Fs^  =  {e*™^m  =  ...,-2,-l,0,l,2,...},  (51) 

and  the  diagonal  matrix  A  defined  by 

A  =  diag{. . . ,  i~^,  1,  i,i^, . . (52) 

Lemma  2.19  Given  A:  >  0,  there  exists  a  function  tfoo  '•  [0, 27r]  C ,  such  that 
for  large  r  >  0,  the  scattered  field  (18)  assumes  the  form 

Proof.  For  arbitrary  integer  m  and  large  z  >  |m|  (see  [3],  pp  364), 

HUz)  =  (•-"  +  0  (;))  ■  (54) 


The  lemma  follows  immediately  from  the  combination  of  (54)  and  (18). 


(46) 

(47) 

(48) 

(49) 
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Definition  2.20  Suppose  that  ^  is  a  scattered  field;  then  the  function 


0oo(^)  =  lim 


2l—  g2(A:r— '7r/4) 
Tzkr  ^ 


is  referred  to  as  the  far  field  ofxf. 

Given  the  incident  angle  fi  G  [0, 27r]  and  the  incident  field 

y)  =  g*^(^cos/3+3/sm/?)  _  ^ikrcos(8-P) 

we  denote  by  ^(r,  d;  fi)  the  corresponding  scattered  field,  and  by 

fi{r,e-,l3) 


r^Qi(kT—KjA) 
V  TTfcr 


(55) 


(56) 


(57) 


the  far  field  of  the  scattered  field  •fi{r,6\fi).  The  function  (57)  is  frequently 
referred  to  as  the  (far-field,  full-aperture)  scattering  amplitude  (see,  for  example, 
[2]  and  [6]).  The  full-aperture  scattering  amplitude  is  related  to  the  scattering 
matrix  (see  Remark  2.5)  via  orthogonal  transforms  specified  in  the  following 
lemma. 


Lemma  2.21  Suppose  that  b  is  a  positive  number  and  S  is  the  scattering  matrix 
corresponding  to  a  scatterer  in  the  disk  D{h).  Suppose  further  that  is 

the  scattering  amplitude.  Then 


m,l 

(58) 

or  equivalently, 

f;oo{OJ)  =  Fi^-A-^-S-A-Fp. 

(59) 

Proof.  It  is  easy  to  verify  that  the  incident  field  (56)  can  be  written 
(17)  with 

in  the  form 

a  =  A-F{fi). 

(60) 

The  corresponding  scattered  field  (18)  therefore  is 

■fi{r,e]fi)  =  Ff^-Ekr-y, 

(61) 

with 

-> 

II 

p 

(62) 

Substituting  (60)  into  (62)  which  is  in  turn  substituted  into  (61),  we 

obtain 

f;{r,e-J)  =  Fi^-Hkr-S-A-F0. 

(63) 
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Now  the  lemma  follows  immediately  from  the  combination  of  (63),  (57)  and  (54). 

□ 

We  say  that  the  scattering  matrix  50,, *  corresponding  to  the  scatterer  q  is 
centered  at  the  origin  (0,0)  since  the  expansions  (17)  and  (18)  of  the  incident 
and  scattered  fields  are  around  the  origin.  Now,  for  the  same  physical  scatterer, 
if  we  shift  the  origin  to  a  point  ^  =  (a,  6),  the  new  scatterer  function  will  be 
q^{x^  y)  =  q{x  —  a,y  —  b)  in  the  new  coordinates.  It  has  compact  support  in  the 
disk  D{A),  A  =  w  +  j^j,  centered  at  the  new  origin.  We  denote  by  5^^^.  the 
scattering  matrix,  corresponding  to  the  same  scatterer,  but  centered  at  The 
following  lemma  is  a  reformulation  of  Lemmas  3.2,  3.3  of  [7]. 

Lemma  2.22  Suppose  that  q  is  the  function  of  a  smooth  scatterer  with  compact 
support  D{w),  that  u  =  (a;i,?/i),  v  =  {x2,y2)  ore  two  points  in  B? ,  and  that 
A  =  zo  +  |u|,  B  =  w  A  \v\.  Suppose  further  that,  corresponding  to  the  same 
scatterer,  S'^k  ond  Sq  ^  are  the  scattering  matrices  centered  at  u  and  v.  Then 

Sl,k  =  T~'^-SXk-T,  (64) 

where  T  :  ^  is  defined  by  the  formula 

J’  _  J?.g*A:[(3/2-yi)-cos(0)-(i2-a:i)-sin(e)].jp-l^ 


where  the  function 


r(«)  =  e 


ik[{y2-yi )  ’  cos(^) - (X2 -XI ) •  sin( ^)] 


(66) 


is  regarded  as  the  diagonal  linear  mapping  F  :  L^[0, 27r]  L^[0,  27r]  defined  by 

the  formula 

(F  .  fm  =  F(^)./(^),  (67) 

for  all  f  e  L^[0,2tt]. 


Remark  2.23  Since  F  and  therefore  T  are  orthogonal  mappings,  two  scattering 
matrices  corresponding  to  the  same  scatterer  but  centered  differently  are  connected 
to  each  other  by  orthogonal  transforms;  therefore,  the  two  scattering  matrices 
contain  the  same  amount  of  information  about  the  scatterer. 


2.5  The  Near  Field 

Given  r  >  0  and  an  incident  field  <f)Q  upon  the  chopped  scatterer  quA)-,  the 
scattered  field  ^  is  smooth  inside  D{r),  continuous  across  the  circle  |x|  =  r,  and 
is  of  the  form  (see  (23)) 

CO 

Y.  (68) 
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outside  D{r).  The  series  (68)  is  absolutely  convergent  for  p  >  r.  We  refer  to 
i^\p=r  3.S  the  near  field.  In  this  subsection,  we  estimate  the  rate  of  convergence  of 
(68)  at  /9  =  r. 

Lemma  2,24  Suppose  that  r  is  a  positive  number  and  that  the  scatterer  q  is 
smooth  on  D{r).  Suppose  further  that  if;  :  E?  C  is  the  scattered  field  corre¬ 
sponding  to  an  incident  field  (f)Q  :  D{r)  — ^  C  upon  the  chopped  scatterer  qo{r)^ 
Then  the  near  field  '0|p=r  has  the  Fourier  series 

oo 

^(r,0)=  ^  (69) 

m=  — OO 


Furthermore,  given  r,  k,  q,  there  exists  c  >  0,  such  that 


brn  ^ 


(70) 


for  all  |m|  >  No{kr). 


Proof.  It  is  well-known  (see  [6])  that  the  scattered  field  is  the  solution  of  the 
Lippmann-Schwinger  equation 

JD{r)  JD{t) 


that  it  is  smooth  on  D{r),  and  that,  given  r,  k,  q,  there  exists  Cj  >  0  such  that 


Halloo  <CiI|^o||2.  (72) 

It  is  also  well-known  that  (see  [1],  [3])  the  free  space  Green’s  function  Gk  can  be 
expressed  as 

i  °° 

G,(x,0  =  --  E  Hm{kp)Jmikg)e^^^^-^\  (73) 

^  m=— CO 

with  X  =  p(cos(0),  sin(0)),  =  ^(cos(t?),  sin(t?)).  Therefore,  the  charge  density 

a  :  D{r)  C  defined  by  the  formula 

a{x)  = -k'^q{x){'il;{x)  + (j)o{x)),  (74) 

is  a  smooth  function,  and  that,  given  r.  A:,  9,  there  exists  C2  >  0  such  that 


Ikiloo  <  C2||<?i'o||2- 
It  follows  from  (71)  and  (74)  that 

V’(aj)  =  /  Gk{x,i)<T{i)di. 

JD{t) 


(75) 


(76) 
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For  X  =  /9(cos(^),sin(0))  and  /?  >  r,  the  substitution  of  (73)  into  (76)  yields 


where 


1  r27r  . 


(78) 


is  the  m-th  Fourier  mode  of  a  on  the  circle  |x|  =  g.  Therefore,  according  to  (75), 
there  exists  C3  >  0  such  that 


|^m(^)|  —  C3||<^o||2? 

for  all  integer  m  and  real  g  E  [0,  r].  Let 

ZTT 

bm  =  ^  /  (l7n{^Q^'Jm{}^Q\HjTi{hT^*Q‘dQ. 

•/  0 

Then 

cz-A\<i>o 


\bm\  <  Jo  \Jmikg)-Hm{kr)\-g-dg. 


(79) 


(80) 


(81) 


For  m  >  No{kr),  according  to  (13)  and  (14),  there  exists  a  positive  number  C4 
such  that  /  \  m 

l^(ie)-J?«(ir)|  <  (^j  .  (82) 

It  follows  immediately  from  (81)  and  (82)  that 

C3-C4-7r||</>o||2  1 


\bm\< 


2  m 


2‘ 


(83) 


Now,  letting  p  r  in  (77)  and  combining  the  result  with  (80)  and  (83),  we 
obtain  (69)  and  (70).  □ 

Remark  2.25  For  \m\  <  z,  it  is  well-known  (see  [3],  pp  366)  that 


(84) 

lu)  ■ 

(85) 

1 

‘(A:r)2/3’ 

(86) 

It  follows  from  (79)  and  (80)  that 

\bm\  <  C3||^o||2 


for  all  |m|  <  kr.  This  means  that  as  a  function  of  m,  the  Fourier  coefficients 
{bm}  of  (69)  remain  fiat  in  magnitude  for  \m\  <  kr.  The  magnitude  begins  to 
decrease  quite  suddenly  when  \m\  becomes  greater  than  kr;  it  decays  to  zero  at 
the  rate  1/m^  for  m  >  No{kr). 
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In  Lemma  2.24,  choose 


(87) 

then  it  follows  from  (22) — (24)  that  the  Fourier  coefficients  {677^,— oo<?7z<oo} 
of  the  near  field  ^|p=r,  ^.re  same  as  the  entries  of  the  the  n-th  column  of  the  matrix 
Hkr'Sr.k-  In  other  words, 

^  i^Hkr ' 771,71  •  (88) 

The  following  lemma  is  a  restatement  of  Lemma  2.24  in  terms  of  the  scattering 
matrix. 

Lemma  2.26  Suppose  that  r  is  a  positive  number  and  that  the  scatterer  q  is 
smooth  on  D[r).  Suppose  further  that  ST,k  is  the  scattering  matrix  corresponding 
to  the  chopped  scatterer  qD{r)-  Then  there  exists  a  positive  number  c  such  that 

\{Hkr-Sr,k)m,n\  <  (89) 

uniformly  for  all  integers  n  and  m  such  that  \m\  >  No{kr). 


Figure  3:  Structures  of  the  Matrices  Sr,k-Hkr,  Hkr-Sr,k  ■ 


The  following  remark  is  a  restatement  of  Remark  2.25  in  terms  of  the  scattering 
matrix. 

Remark  2.27  When  m  becomes  greater  in  absolute  value  than  kr,  the  elements 
on  the  m-th  row  of  Hkr'Sr,k  begin  to  decrease;  they  decay  uniformly  at  the  rate 
1/ivf  for  m  >  No{kr).  If  we  omit  the  small  entries  in  the  part  of  the  matrix 
Hkr'Sr,k  whose  row  index  m  is  greater  than  No{kr)  in  absolute  value,  it  follows 
from  Remark  2.14  Hkr'Sr,k  is  approximately  a  square  matrix  of  dimension 
2-A^o(^’')  (see  Figure  3).  It  follows  immediately  from  Lemmas  2.17  and  2.18  that 

{Sr^k-HkrY  =  T-Hkr-Sr,k-  (90) 

Therefore,  {SrxHkrY  is  approximately  a  square  matrix  of  dimension  2-No{kr), 
aUid  so  zs  S .^^k* H kv  * 
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Figure  3  shows  the  structures  of  the  two  matrices  Sr,k-Hkr  and  Hkr-Sr,k',  compare 
with  Figure  2.  Only  the  first  quadrant  of  the  matrices  (entries  with  both  row  and 
column  indices  positive)  is  depicted  here.  In  the  figure,  the  section  of  a  matrix 
labeled  zero  have  entries  small  compared  to  entries  in  the  square  submatrix  of 
dimension  2-No{kr). 


2.6  Continuation  Method  for  Nonlinear  Problem 


For  simplicity,  we  will  restrict  our  discussion  to  the  case  of  finite  dimensional 
space.  Let  us  consider  a  nonlinear  mapping  P  :  i?"  x  [0, 1]  P"  and  the 
solution  of  the  problem 

P(x,l)  =  0.  (91) 

Suppose  that  for  every  A  G  [0, 1],  there  exists  a  unique  solution  xx  G  P”  to  the 
problem 

P(a:A,A)  =  0.  (92) 

Suppose  further  that  xx  depends  smoothly  on  A,  and  that  the  mapping  P  is 
smooth.  Finally,  suppose  that  xq  is  known.  Then  there  is  a  procedure  referred 
to  as  the  continuation  method  which  obtains  the  solution  xi  using  the  “initial 
solution”  xo,  by  recursively  solving  a  series  of  linearized  problems.  This  simple 
scheme  is  as  follows. 

Suppose  that  at  some  A  G  [0,1)  we  have  obtained  the  solution  xx-  For  a 
sufficiently  small  h  >  0,  we  intend  to  obtain  xx+h,  the  solution  of  the  equation 

P{xx+k,X  +  h)  =  0.  (93) 


Subtracting  (92)  from  (93),  we  have 

idP{xx,\) 


V 


(xa+a  -  Xx)  =  +  0{h^). 


dx  j  .  dX 

In  other  words,  up  to  the  second  order  of  h,  the  perturbation 

6x  =  xx+h  —  Xx 


(94) 


(95) 


is  a  solution  of  a  linear  problem.  If  we  further  assume  that  the  n  x  n  matrix 

_ ^ 

dx 

is  not  singular,  the  increment  8x  can  be  determined  up  to  second  order  of  h. 

Remark  2.28  In  most  applications,  once  the  second  order  approximation  of 
Xx+h  is  obtained,  no  further  attempt  is  made  to  compute  the  exact  xx+k-  in¬ 
stead,  the  recursion  in  X  proceeds  to  compute  the  next  solution  xx+2h  from  the 
approximate  solution  Xa+a  j'^st  obtained. 
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Obviously,  the  linearization  procedure  is  a  standard  perturbation  analysis  on 
the  parameter  A.  Sometimes,  it  is  more  convenient  and  advantageous,  both 
computationally  and  analytically,  to  carry  out  the  perturbation  analysis  in  a 
slightly  different  way.  Suppose  we  do  not  wish  to  evaluate  the  term 


gP(xA,  A) 
dX 

on  the  right  hand  side  of  (94).  Given  xx,  let  us  compute 

c  =  P{x\,  X  h). 

Subtracting  (98)  form  (93),  we  obtain 

(dP{xx,X  +  h) 


V 


dx 


■Sx  —  — e  +  0{h^). 


(97) 


(98) 


(99) 


It  is  easy  to  verify  that  e  is  of  order  h.  As  before,  a  second  order  approximation 
of  the  solution  xx+h  can  be  obtained  from  xx  by  solving  a  linear  problem. 


2.7  The  Solution  of  a  Linear  ODE  of  Matrix 

Suppose  that  A(r),  B{r)  and  C{r)  are  three  n  x  n  matrices  depending  contin¬ 
uously  on  r  G  [0,1].  Let  us  consider  an  ordinary  differential  equation  of  the 
form 


5'(r)  =  A(r)-S'(r)  +  5(r)-5(r)  +  C{r).  (100) 

Now,  we  wish  to  express  in  close  form  the  solution  S{r)  at  an  arbitrary  r  €  (0, 1] 
for  prescribed  initial  value  5(0). 

Lemma  2.29  Suppose  that  the  nxn  matrix  P{r)  is  continuous  for  all  r  €  [0, 1]. 
Suppose  further  that  a  <  b  are  two  real  numbers  in  [0, 1].  Finally  suppose  that  m 


is  a  positive  integer,  and  that  h  =  {b  —  a) fm.  Then  the  two  limits 

EL{a,h-,P)  =  +  h-P(b)){I  +  h-P{b  —  h))  X  . . .  X 

iI  +  h-Pia  +  h)){I+h-P{a)),  (101) 

ER{a,h-,P)  =  +  h-P{a)){I  +  h-P{a  +  h))  X  . . .  X 

(I +  h-P{ib-h)){I  +  h-Pib))  (102) 

exist.  Moreover,  for  arbitrary  real  numbers  rj,  r2,  ra  G  [0, 1], 

EL{ri,rz]P)  =  EL{r2,r3]P)EL{ri,r2]P),  (103) 

pR{ri,r3-,P)  =  £R(ri,r2;P)Efl(r2,r3;P).  (104) 
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Finally^ 


"  f;i(o.6;/>)-P(o),  (105) 

=  P(h)-Ei,(a,h-,P),  (106) 

=  P{a)-ER(i,kP),  (107) 

=  ER(a,b-,P)-P(b).  (108) 

The  proof  of  Lemma  2.29  is  trivial,  and  is  omitted  here. 


Remark  2.30  If,  for  any  real  numbers  c,  d  ^  [0, 1],  P{c)  commutes  with  P{d), 
then 

Elia,  fe;  P)  =  ER{a,  6;  P)  =  exp 

Remark  2.31  Under  the  conditions  of  the  preceding  lemma,  it  can  be  easily 
shown  that 

£i(a,i;P)  =  {l-^P(l>)j  (/  +  5-P(<'-A))  X 

I-^Pib-h))  (^/ +  |p(6  -  2/.)^  X  . . .  X 

-  jP(a  +  i)  j  ^/ +  |p(a)  j  +  0(5”),  (110) 

and 

En(a,h-,P)  =  (l+  ^P(a)'j  f  /  -  |p(a  +  *)^  x 

/+|p(<.  +  5)jf/- jP(o  +  2A)^  X...X 

(l+^P{b-h)'j(^l-^P{b)'j  +0(5”).  (Ill) 

The  following  lemma  is  an  immediate  consequence  of  the  preceding  one. 


Lemma  2.32  Suppose  that  24(r),  B{r)  and  C{r)  are  three  n  x  n  matrices  de¬ 
pending  continuously  on  r  £  [0, 1] .  Then 

S{r)=  r  EL{T,r;AyC{T)-ER{r,r-,BydT  (112) 

Jo 
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is  the  solution  of  the  initial  value  problem 


(113) 

(114) 


5'(r)  =  A{ryS{r)  +  Sir)-B{r)  +  C{r), 

5(0)  =  0. 

3  Heisenberg’s  Uncertainty  Principle  and  Re¬ 
cursive  Linearization 

It  turns  out  that  the  ill-posedness  of  the  inverse  scattering  problem  can  be  ben¬ 
eficially  used  to  solve  it.  It  means  that,  due  to  ill-posedness  of  the  problem,  not 
all  equations  in  the  nonlinear  system  (see  Remark  2.15)  are  strongly  nonlinear, 
and  that  when  solved  recursively  in  a  proper  order,  they  can  be  reduced  to  a 
collection  of  linear  problems.  In  this  section,  we  reformulate  the  ill-posedness 
and  the  inverse  scattering  problem,  and  present  an  inversion  algorithm.  More 
specifically,  in  Section  3.1,  we  examine  and  reformulate  the  ill-posedness  of  the 
inverse  scattering  problem  in  the  special  case  of  weak  scattering.  In  Section  3.2, 
we  briefly  and  informally  describe  the  recursive  linearization  method.  In  Section 
3.3,  we  present  the  Heisenberg’s  Uncertainty  Principle  for  the  scattering  problem, 
which  we  shall  use  in  Section  3.4  to  reformulate  the  inverse  problem.  The  details 
of  the  inversion  algorithm  will  be  described  in  Section  3.6. 

3.1  A  Special  Case 

When  5  or  A:  is  small,  the  scattered  field  is  weak,  and  the  inverse  scattering  prob¬ 
lem  becomes  essentially  linear.  In  this  subsection,  we  examine  this  special  case 
and  make  necessary  connections  to  the  general  case  where  the  inverse  problem  is 
nonlinear. 

Nowhere  does  the  ill-posedness  of  the  inverse  scattering  problem  become  more 
manifest  than  in  the  case  of  weak  scattering.  As  is  well-known  (see,  for  example, 
[8],  [12]),  when  ^  <C  1  or  A:  1,  the  scattered  field  is  weak,  where  the  Born 

Approximation  to  the  far  field  ipooiG,0)  (see  Section  2.4)  can  be  written  in  the 
form 

with  an  error  of  0{q^)  if  q  is  small,  or  of  0(A;'*  log^(A:))  if  k  is  small.  In  other 
words,  under  the  assumption  of  weak  scattering,  the  far  field  ■0oo(^,/5)  depends 
on  q  essentially  linearly,  and,  up  to  a  higher  order  error  and  a  scaling,  is  the 
Fourier  transform  of  q 

qm,n  =  y)e^^^^^+^^Uxdy  (116) 
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(117) 

(118) 


with  the  pair  of  real  numbers  m  and  n  given  by 

m  =  A:(cos/3  —  cos^), 
n  =  k{sin  13  —  sinO). 


Therefore,  the  full-aperture  far-field  measurements 

{  tor  all  (ej)  e  10,2,r]  X  (0,27r]  } 


(119) 


are  the  Fourier  transform  qm,n  for  those  points  (m,  n)  filling  the  entire  disk  D{2k), 
which  we  refer  to  as  the  Fourier  aperture  of  radius  2k  (see  Definition  2.4  and 
Figure  5).  With  such  measurements,  the  scatterer  q  can  be  determined,  obviously, 
with  the  resolution 


27r 


a  = 


radius  of  Fourier  aperture 


—  -lx 

2k  ~  2^’ 


(120) 


where 


is  the  wavelength.  We  consequently  have 


(121) 


Lemma  3.1  (Uncertainty  Principle,  Small  q)  Suppose  that  q  is  small.  Then 
from  the  far-field  measurements,  we  cannot  determine  features  of  the  scatterer 
that  are  less  than  half  a  wavelength. 


Remark  3.2  Lemma  3.1  is  a  reformulation  of  the  ill-posedness  of  the  (linear) 
inverse  problem.  It  explicitly  specifies  the  null  space  of  the  linear  operator  (115) 
which  maps  the  scatterer  q  to  the  scattering  measurements  V’oo-  the  Fourier  modes 
of  q  higher  than  2k  are  not  observable  in  the  measurements,  and  thus  cannot  be 
determined. 


Remark  3.3  Formulae  (115)— (121)  are  also  valid  for  small  k;  therefore.  Lemma 
3.1  implies  that  at  a  sufficiently  low  frequency  only  ^(0,0)  (the  average  of  the 
scatterer)  can  be  determined  from  the  measurements. 


The  scattering  matrix  Sr,k  is  small  if  scattering  is  weak.  Thus,  for  small  q  or 
small  k,  the  Riccati  equation  (34)  can  be  linearized  by  omitting  terms  from  the 
right  hand  side  of  (34)  quadratic  and  cubic  in  q  and  Sr,k-  The  solution  of  the 
linearized  equation 

ds^  = 


dr 


is  given  by  the  formula 


^■u^,k 


ixk^ 


I  Jkr'^r'Jki 

Jo 


•r-dr. 


(122) 


(123) 
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Since  is  an  approximation  to  the  scattering  matrix  S^^k,  it  should  be  con¬ 
nected  to  the  Born  Approximation  (115)  due  to  formula  (58).  The  following  is  a 
restatement  of  Lemmas  2.18  of  [11]  in  terms  of  the  scattering  matrix. 

Lemma  3.4  Suppose  that  is  the  solution  of  (122),  Sr,k  is  the  solution  of  the 
Riccati  equation  (Sf).  Then  there  exists  a  constant  c  >  0,  such  that  for  k  >  0 
and  all  r  G  [0,  z:u]j 

-  s.®  lb  <  C  •  (i  ■  |log(m  ■  ||?|U)' .  (124) 

Moreover, 

J  Q 

3.2  The  Inversion  Algorithm,  An  Informal  Description 

Let  us  denote  by 

P{q,  k)  =  S^,k  (126) 

the  system  of  nonlinear  equations  for  the  inverse  scattering  problem.  In  this 
subsection,  we  briefly  describe  a  simple  procedure  that  solves  the  inverse  problem. 

For  a  given  precision  e  >  0  and  frequency  k  >  0,  there  should  be  infinite 
number  of  forward  models  q  that  satisfies  (126)  to  the  prescribed  precision,  due 
to  the  ill-posedness  of  the  problem.  We  choose  from  them  the  most  smooth  one 
and  denote  it  by  qk-  Therefore,  to  the  given  precision,  the  inverse  problem  can 
be  reformulated  as 

P{qk,  k)  =  (127) 

We  expect  that  qk  lives  in  a  finite  dimensional  subspace  of  just  as  it  is  the 
case  when  q  is  sufficiently  small;  the  Fourier  transform  of  q^  is  essentially  zero 
outside  the  disk  D(2k)  (see  Remark  3.2).  The  inversion  algorithm  is  a  recursive 
linearization  procedure  which  recovers  qk  from  small  k  to  large  k. 

For  sufficiently  small  k,  according  to  Remark  3.3,  qk  lives  in  a  one-dimensional 
subspace;  it  is  the  average  of  q.  Moreover,  the  equation  (127)  is  linear  to  the 
prescribed  precision  e,  and  therefore  can  be  solved  in  the  least-squares  sense  to 
obtain  qk- 

Now  suppose  that  qk  depends  continuously  on  k,  and  that  we  have  recovered 
qk  at  some  A:  >  0.  Then  the  standard  procedure  of  continuation  (see  Section  2.6) 
can  be  used  to  recover  qk+sk  by  solving  a  linear  problem  for  the  perturbation 

^9  =  qk+sk  -  qk-  (128) 

Consequently,  the  inverse  scattering  problem  (126)  can  be  solved  up  to  any  given 
frequency  k  and  to  the  prescribed  precision  e  >  0,  provided  that  the  scattering 
data  (20)  are  available. 
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In  the  next  several  subsections,  we  characterize  more  precisely  the  finite  di¬ 
mensional  subspace  in  which  qu  resides,  and  describe  the  inversion  algorithm  in 
detail. 


3.3  Uncertainty  Principle,  the  General  Case 

According  to  Lemma  3.4,  when  q  is  small  and  up  to  the  second  order  of  the  small¬ 
ness,  the  knowledge  of  the  scattering  data  is  equivalent  to  the  knowledge  of 
the  Fourier  modes  of  q  in  the  aperture  D{2k).  It  turns  out  that  when  q  is  not 
small,  the  above  statement  still  essentially  holds.  In  other  words,  S^,k  contains 
information  of  the  Fourier  modes  of  q  essentially  in  D{2k).  In  this  subsection, 
we  wish  to  make  this  assertion  more  precise,  and  to  show  it  is  indeed  the  case. 

In  the  evaluation  of  the  right  hand  side  of  the  Riccati  equation  (34),  Jkr  + 
Sr,k-Hkr,  Hkr-Sr,k  +  Jkr,  according  to  Remarks  2.13  and  2.27,  essentially  are  two 
square  matrices  of  size  2‘No{kr).  Therefore,  the  operation 

{Jkr  “1"  Sr,k‘Hkr)'qr'{Hkr'ST,k  "h  Jkr} 


on  qr  (see  (33))  is  a  process  of  low-pass  filtering  on  the  scatterer  q,  as  depicted 
in  Figure  4.  At  a  given  frequency  k  and  on  the  circle  |a;|  =  r,  only  the  Fourier 
modes 

{  q,^{r)  I  |m|  <  2-No{kr)  }  (130) 

essentially  participate  in  the  operation;  higher-frequency  modes  of  the  scatterer 
are  filtered  out  in  the  process.  The  relatively  low-frequency  angular  Fourier 
coefficients  (130)  at  r  are  therefore  picked  up  in  the  integration 


Sr.,k  =  ^  r(Jkr  +  Sr,k-Hkr)-qriHkr-Sr,k  +  JkrYr-dr,  (131) 
2  Jo 

and  encoded  in  the  scattering  data  We  thus  conclude  that  the  scattering 

data  contain  insignificant  information  of  the  higher-frequency  angular  Fourier 
coefficients 

{  qm{r)  1  |m|  >  2-No{kr)  }  (132) 

of  the  scatterer  for  all  r  €  [0, 137].  In  other  words,  at  each  r  >  0,  the  resolution  of 
the  scatterer  on  the  circle  |a;]  =  r  is 


27rr  ^  ^ 
2No{kr)  k  2 


(133) 


That  is,  the  angular  resolution  is  about  half  a  wavelength.  Features  of  q  smaller 
than  half  a  wavelength  in  the  angular  direction  contribute  considerably  weakly  to 
the  scattering  data;  the  smaller  the  features  become,  the  weaker  the  contributions 
are,  and  the  more  difficult  it  becomes  to  recover  these  small  features.  The  above 
discussion  can  be  summarized  in  the  following  lemma. 
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Figure  4:  Process  of  Low-pass  Filtering  of  the  Scatterer. 


Lemma  3.5  (Uncertainty  Principle,  Angular  Direction)  Suppose  that  the  smooth 
scatterer  q  has  compact  support  D{w),  and  that  is  the  corresponding  scatter¬ 
ing  matrix.  Then  it  is  increasingly  difficult  to  determine  from  Szp,k  small  features 
of  q  in  the  angular  direction  as  their  sizes  become  increasingly  less  than  half  a 
wavelength. 

Now,  we  show  that  the  resolution  of  the  scatterer  in  an  arbitrary  direction  is  also 
about  half  a  wavelength.  Given  f  =  {a,b)  €  R^,  let  A  =  w  +  |^|,  and  consider 
the  scattering  matrix  S^  j^,  centered  at  ^  and  corresponding  to  the  scatterer  q 
(see  Section  2.4).  According  to  Remark  2.23,  and  the  scattering  data  S^^k 

contain  the  same  amount  of  information  about  the  scatterer;  therefore  f.  can 
be  used  as  the  scattering  data  to  recover  the  same  scatterer,  now  represented  by 
the  function  q^{x,y)  =  q(x  —  a,y  —  b).  According  to  Lemma  3.5,  the  angular 
resolution  provided  by  the  scattering  data  is  about  half  a  wavelength.  In 
other  words,  the  resolution  of  q^  on  any  circle  centered  at  f  is  about  half  a 
wavelength.  Therefore,  for  a  given  point  x  inside  the  scatterer  and  for  a  given 
direction  r,  we  can  always  choose  the  center  which  is  sufficiently  far  from  the 
location  of  the  scatterer,  such  that  there  is  one  circle  centered  at  ^  which  passes 
through  X  in  direction  r.  Since  the  resolution  of  the  scatterer  on  this  circle  is 
about  half  a  wavelength,  and  since  contains  the  same  amount  of  information 
about  q  as  the  original  scattering  data  Fc?,*  does,  we  obtain: 

Lemma  3.6  (Uncertainty  Principle)  Suppose  that  the  smooth  scatterer  q  has 
compact  support  Dffi),  and  that  S.^^k  is  the  corresponding  scattering  matrix. 
Then  it  is  increasingly  difficult  to  determine  from  S^^k  small  features  of  q  as  their 
sizes  become  increasingly  less  than  half  a  wavelength.  In  other  words,  contained  in 
the  measurements  S^^k  are  essentially  the  Fourier  modes  of  q  inside  the  Fourier 
aperture  D{No{2k)). 
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Denote  by  ^(2^)"^  the  Fourier  aperture  where  the  Fourier  modes  of  q  can  be  re¬ 
covered  from  the  measurements  in  a  well- conditioned  procedure.  Presently, 
no  detailed  characterization  of  D{2k)'^  is  available  except  that  it  belongs  to 


D{No{2k)). 


Figure  5:  Fourier  apertures  D(2k)  and  D{2k)'^ . 

Remark  3.7  The  Uncertainty  Principle  is  an  equivalent  formulation  of  the  ill- 
posedness  of  the  inverse  scattering  problem:  small  features  of  the  scatterer  belong 
to  the  virtually-null  space  (see  Remark  2.12)  of  the  nonlinear  mapping  of  the 
inverse  scattering  problem;  they  are  essentially  not  observable  in  a  scattering 
experiment. 


3.4  Reformulating  Scattering  Problem 

Denote  by  qk  the  low-frequency  part  of  q  in  the  Fourier  aperture  D{2k)'^ .,  so  that 


J  9m, n,  (m,n)  €  jD(2A;)+, 

\  0,  (m,n)  ^  D{2kY . 


(134) 


The  goal  of  inversion,  in  the  lights  of  Lemma  3.6,  is  to  stably  obtain  qk  within 
a  reasonable  precision.  By  definition,  qk  is  the  only  component  of  q  that  is 
observable  in  the  scattering  data;  consequently  the  original  forward  scattering 
model  q  can  be  replace  by  qk  without  (essentially)  changing  the  measurements. 
We  therefore  can  reformulate  the  equation  (34)  as 


^  =  '-^k%Jkr  +  Sr,k-Hkr)iqk)AHkr-Sr,k  +  Jkr).  (135) 

dr  I 

Definition  3.8  To  a  scattering  experiment  at  frequency  k,  a  scatterer  q  is  said 
to  look  (essentially)  the  same  as  a  scatterer  q  if  they  produce  essentially  the  same 
scattering  measurements  in  the  experiment. 
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Definition  3.9  A  forward  model  q  is  said  to  be  observable,  or  an  observable  part 
of  the  original  scatterer  q,  to  a  scattering  experiment  at  frequency  k,  if  it  looks 
the  same  as  the  original  q,  and  its  Lf  norm  is  the  least  among  those  that  look  the 
same  as  q. 

Remark  3.10  At  frequency  k,  q^  looks  the  same  as  the  original  q  to  a  full- 
aperture  experiment,  it  is  also  observable  to  the  full-aperture  experiment.  On  the 
other  hand,  in  an  experiment  of  limited  aperture,  qk  may  not  be  the  observable 
forward  model,  but  it  looks  the  same  as  the  observable. 

3.5  Continuity  of  in  Frequency  k 

In  this  subsection,  we  argue  that  qk,  the  observable  part  of  q  at  frequency  k, 
depends  on  k  continuously.  This  is  certainly  true  in  the  special  case  of  small  q. 
There,  the  observable  part  of  q  is  the  Fourier  modes  of  q  in  aperture  D{2k)  (see 
Section  3.1).  Therefore,  new  Fourier  modes  added  to  qk+5k  are  those  qm.,n  in  the 
annulus 

A{k,8k)  =  {  (m,n),  2k  <  Vm^  n'^  <  2{k  +  Sk)  }.  (136) 

Consequently,  the  perturbation  in  qk,  due  to  the  small  perturbation' in  k,  is  small; 

\\qk+6k  -  qkh  =  Ikfc+ifc  -  qfh  =  /  gm,n  ■  dm  -  dn  =  0{8k).  (137) 

JA{k,6k) 

Assuming  the  well-posedness  of  the  initial  value  problem  (see  (34),  (26))  of 
the  Riccati  equation,  we  further  argue  that  the  dependence  of  qk  on  k  is  also 
continuous  in  the  general  case  where  q  is  not  small.  This  well-posedness  means, 
in  particular,  that  the  scattering  matrix  S'r,*  is  a  smooth  function  of  k,  that  its 
virtual  rank  and  structure  of  essentially  nonzero  entries  (see  Figure  1)  depend  on 
k  smoothly.  Therefore  the  amount  of  information  the  process  (129)  acquires  from 
q  depends  on  k  smoothly.  We  summarize  the  above  discussion  as  an  Observation 
for  later  reference. 

Observation  3.11  To  a  full-aperture  experiment,  the  observable  scattering  model 
qk  depends  continuously  on  k  in  the  Lf  norm. 

We  wish  to  carry  this  point  further  to  the  case  of  limited- aperture  measurements. 
Denote  by  qk,i  the  observable  part  of  q  corresponding  to  an  experiment  of  a  limited 
aperture.  Usually,  qk,i  is  not  the  same  as  qk,  and  therefore,  due  to  Definition  3.9, 

lkfc,i||2  <  hkh-  (138) 

We  postulate  that  Observation  3.11  is  also  valid  for  scattering  experiments  with 
limited  aperture.  This  has  been  observed  in  our  numerical  experiments,  and  can 
be  proved,  again,  in  the  special  case  of  small  q. 
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3.6  A  Recursive  Linearization  Algorithm 

Suppose  that  a  set  of  full-aperture,  full- bandwidth  scattering  data  (see  (20))  are 
given,  we  present  in  this  subsection  a  stable  method  for  the  solution  of  the  inverse 
scattering  problem.  There  are  two  approaches  to  the  description  of  the  method: 
one  is  based  on  the  Lippmann- Schwinger  equation  (see  [12]),  the  other  on  the 
Riccati  equation,  which  has  been  numerically  implemented  (see  Section  4),  and 
which  we  wish  to  present  here. 

Let  us  again  consider  the  nonlinear  mapping  (see  Remark  2.15)  which  maps 
the  scatterer  q  to  the  scattering  data  (see  (20)) 

{  S^,k  1  0  <  A;  <  oo  },  (139) 

defined  by  the  initial  value  problem  (see  Section  3.4  and  (26)) 

S'r,k  =  -Yk\Jkr  +  SrxHkr)iqk)AHkr-Sr,k  +  Jkr),  (140) 

So,k  =  0.  (141) 

Discretizing  the  /c-variable  with  nodes  h,  kz,  ...  (see  Figure  6),  we  now 
describe  a  procedure  which  recursively  determines  qk  sA,  k  =  kj  for  j  =  1,2,... 
Indeed,  for  sufficiently  small  fei,  the  relationship  between  qk^  and  Sr,ki  becomes 

0  k 

I - 1 - 1 - i - i - 1 - 1 - 1 - 1 - 1 - 1 - i - i - i - 1 - i - 1 - ► 

h  k2  h 

Figure  6:  Computational  Grid  in  the  Frequency  Space. 

essentially  linear  (see  Section  3.1),  and  the  problem  (140),  (141)  can  be  replaced 
by  the  linear  one  (see  Lemma  3.4): 

SU.  =  (142) 

So.k,  =  0  (143) 

The  solution  to  this  linear  initial  value  problem  is  obviously  given  by  the  formula 

=  V/o  •’•^AlbA-Ar-r-dr.  (144) 

Now,  for  the  given  scattering  data  S.co,ki,  the  scatterer  qk-y  can  be  obtained  by 
solving  the  linear  problem  (144). 

Remark  3.12  In  practice,  only  an  approximate  qk-,  is  needed  to  start  the  recur¬ 
sive  procedure,  namely,  to  determine  qk^-  Therefore,  ki  can  usually  be  chosen 
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quite  large,  and  consequently,  the  linear  equation  (142)  is  not  a  fine  approxima¬ 
tion  to  the  Riccati  equation  (I4O).  We  find  in  numerical  experiments  that  the 
lowest  frequency  ki  can  be  chosen  such  that  its  corresponding  wavelength  is  about 
the  size  of  the  scatterer. 

By  induction,  suppose  that  the  scatter  qj^  have  heen  recovered  at  some  k  >  0,  and 
0  k  k 

I - i - i - 1 - 1 - 1 - 1 - i - 0 - 1 - i - ^ - 1 - 1 - 1 - 1 - 1 - ► 

Figure  7:  Update  from  k  to  k. 

that  A;  >  0  is  slightly  greater  than  k.  We  intend  to  determine  q^,  or  equivalently, 
to  determine  the  perturbation 


Sq=^  qk-  q-k-  (145) 

This  can  be  achieved  by  employing  the  perturbation  analysis  of  Section  2.6  where 
the  continuation  parameter  A  is  now  the  frequency  k.  Following  the  procedures 
described  in  (98)  and  (99),  we  solve  at  the  frequency  k  the  forward  scattering 
problem 

^r,k  ~  ~^k^{JkT  +  Sr,k'^kr)'{qk),.'{Hkr'Sr,k  +  Jkr)^  (146) 

So,k  =  0,  (147) 

corresponding  to  the  scatterer  q^.  As  a  result  of  the  forward  solve,  we  obtain 
ST,k  for  all  r  €  [0,  ct],  which  will  be  used  later  in  the  linearized  equations  for  the 
perturbation  6q  (see  equations  (151),  (153)). 

Remark  3.13  On  the  assumption  that  the  initial  value  problem  of  the  Riccati 
equation  (see  (26),  (34))  is  well-posed,  we  observe  that  the  scattering  matrix  Sr ^k} 
the  solution  of  (146)  and  (14^)  with  the  scatterer  function  q-^,  is  different  from 
but  close  to  Sr,k;  for  the  latter  is  the  solution  of  the  same  equations 

^r,k  —  ~^^^{JkT  +  Sr,k-Hkr)‘iqk)r'{Hkr-Sr,k  +  Jkr),  (148) 

=  0,  (149) 

with  a  different  scatterer  function  qk  close  to  qj.. 

Subtracting  (146)  from  (148),  and  omitting  the  second  order  smallness  in  6q  and 
in 

{,6Sfi^k  ~  Sr,k  Sr^kj  (150) 
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we  obtain  the  linear  ordinary  differential  equation 

+  {^S)r,k'Hkr-i<lk)r'{Hkr-Sr,k  +  Jkr)  + 

+  {Jkr  +  Sr,k-Hkr)-{Sq)r‘{Hkr-ST,k  +  Jkr)}  (151) 

for  6S,  with  its  initial  value 

(55)o,fc  =  0.  (152) 

According  to  Lemma  2.32,  the  linear  initial  value  problem  (151),  (152)  has  the 
formal  solution 

(55),., fc  =  /  EL{T,r]Py)-{Jkr  +  Sr,k-HkT)  ^ 

Z  Jo 

(M)^iHkr-$r,k  +  JkryERir,  r;  P^yr-dr  (153) 

where  Pa,(r)  :  Xkr  ^  Xkr,  Py{r)  :  Ykr  Ykr  are  defined  by  the  formulae 

p,{r)  =  jPrHkr-(%}AHkr-Sr,k  +  Jkry  (154) 

Pyir)  =  jk^iJkr+SrrHkryQy-Hkr.  (155) 

In  particular,  &t  r  =  w,  (153)  becomes  a  system  of  linear  equations  for  6q  : 

r  EL{r,W,Pyy{Jkr  +  Sr, k-Hkry(^l  X 
Z  Jo 

{Hkr-Sr,k  +  JkryEAv,  W,  P^T-dr  =  {8S)^,k,  (156) 

with  the  right  hand  sides  (55)c7,fc  given,  and  the  coefficients  El,Er,S  known. 
Denote  by  L{Xkw  ifc®)  Ibe  linear  space  of  all  linear  mappings  from  Xko  to 
Ykm,  and  by  Ck  '■  L^[D(r<7)]  L{Xkv,  ^  YkA  the  linear  operator  defined  by 
(156).  The  linear  equation  (156)  can  be  rewritten  as 

Ck{Sq)  =  (55),., fc.  (157) 

The  linear  equations  can  be  solved  (see  Remark  3.14)  for  Sq,  and  the  scatterer 

qk  can  be  obtained  from  the  previously  recovered  scatterer  qj.  via  (145). 

Remark  3.14  The  virtual  rank  (see  Remark  2.12)  of  the  linear  operator  Ck  is 
finite  due  to  the  ill-posedness  of  the  inverse  scattering  problem  (see  Lemma  3.6 
and  Remark  3.7).  Therefore,  (157)  is  solved  as  a  least-squares  problem  of  finite 
dimensions  to  yield  the  solution  6q.  Our  numerical  experiments  show  that  in 
fact  only  an  roughly  approximate  qk  is  needed  to  continue  the  up  recursion  in  k, 
namely,  to  stably  determine  the  scatterer  at  the  next  higher  frequency. 
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4  Implementations  of  the  Recursive  Procedure 

In  this  section,  we  discuss  the  discretization  of  the  spatial  variables  (r,  0),  the 
treatment  of  the  scattering  matrix  S,  the  numerical  computation  of  the  Riccati 
equation  (146),  the  evaluation  of  the  linear  operators  El,  Er,  the  formation  of 
the  linear  system  (156),  and  the  least-squares  solution  of  it. 


4.1  Discretizing  the  Independent  Variables  {r,9) 

4.1.1  Discretizing  the  Azimuth  9 

Given  a  radius  r  >  0  and  an  even  number  A  >  0,  we  denote  by 

{  I  ^  }  (15S) 

the  equispaced  points  on  the  circle  {  {r,9)  \  6  e  [0, 27r]  }  (the  value  of  N  will  be 
specified  later  in  this  subsection),  so  that  functions  on  the  circle  are  represented 
by  their  values  at  these  points.  In  particular,  the  scatterer  q  and  its  perturbation 
8q  on  the  circle  are  understood  as  real  valued  vectors  of  dimension  iV;  the  linear 
diagonal  operator  qr  (see  Remark  2.8)  is  now  a  diagonal  matrix  of  dimension 
N\  the  Fourier  transform  E  and  its  inverse  F~'^  are  understood  as  the  discrete 
Fourier  transforms  (DFT)  of  dimension  A;  and  the  linear  operator 

qr  =  F-qr-F-^  (159) 

(see  (32))  is  regarded  as  a  matrix  of  dimension  N.  A  sequence  a  =  {om}  G  Xkr 
is  truncated  and  rearranged  in  the  DFT  order 

{  ao,  Oi, . . . ,  a;v/2,  Q;_;v/2-i-i,  •  •  •  lOf-i  };  (160) 

so  truncated  and  rearranged  are  the  vectors  in  Ykr,  the  matrices  Jkr,  Hkr  and 
Sr,k-  We  will  refer  to  the  central  entries  of  the  vector  (160)  as  the  high-frequency 
entries.  The  high-frequency  entries  of  the  scattering  matrix  are  those  in  the 
center  rows  and  columns.  A  vector 

{  /^O,  /^l,  •  •  •  ,  9m/2,  /^-M/2+l5  •  •  •  5  /5_i  }  (161) 

of  dimension  M  <  N  {M  even)  can  be  added  to  a  by  firstly  zero  padding  j3 

{  /3o,/?i,-..,/?m/2,0,  ...0,^_m/2+i,---,/5-i  },  (162) 

and  then  carrying  out  the  regular  addition  of  two  vectors  of  the  same  size.  Sim¬ 
ilarly,  a  scattering  matrix  of  dimension  M  can  be  viewed  as  of  dimension  N  by 
zero  padding. 
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Remark  4.1  When  matrices  of  different  dimensions  appear  in  an  arithmetic 
operation  (see,  for  example,  (186)),  the  smaller  matrices  are  first  zero  padded  to 
the  maximum  dimension;  the  final  result  is  a  matrix  of  the  maximum  dimension. 

In  the  rest  of  this  paper,  we  denote  by  ne(r)  =  N .  Numerical  experiments 
show  that  ne{r)  can  be  chosen  between  2rk  and  Zrk,  which  is  to  say,  2  to  3 
points  per  wavelength  along  the  arclength  of  the  circle  |x|  =  r.  In  our  numerical 
experiments, 

n0{r)  >  2.8rk.  (163) 

4.1.2  Discretizing  the  Radius  r 

Over  the  interval  [0,ru],  we  employ  two  sets  of  equispaced  computational  grids: 

{rj=j-K,  j  =  0,1,..., Ur,  hr  =  —  },  (164) 

Ur 

{pm=^m-hp,  m  =  0,1,..., Up,  hp  =  —  },  (165) 

Tip 

with  integers  >  Up.  The  first  set  of  grid  is  used  for  the  solution  of  initial  value 
problem  (146)  and  (147).  The  second  set  is  for  discretization  of  the  integral  in 
equation  (156).  Our  experiments  show  that,  with  the  second  order  ODE  solver 
of  Section  4.2, 

Ur  ~  10^^,  (166) 

TT 

namely,  the  grid  is  about  ten  points  per  wavelength  over  the  interval  [0,0?].  For 
convenience  in  computation,  the  ratio  n^/up  is  chosen  as  an  integer,  so  that 

{  Pm  }  C  {  r,-  }.  (167) 

The  ratio  is  2  or  3  in  our  numerical  experiments. 

4.2  Solving  the  Forward  Scattering  Problem 

The  initial  value  problem  (146)  and  (147)  of  the  Riccati  equation  is  solved  by  an 
second  order,  implicit,  alternating  scheme.  Assuming  that  j  is  the  step  counter, 
we  initially  set  ^’  =  0,  r  =  0,  and 

So,k  =  0.  (168) 

For  j  =  1,2, ...  ,nr,  S  is  updated  from  to  rj  by  the  following  procedures 
(with  r  =  [rj  +  rj_i)/2):  if  j  is  odd, 

Z'lT'T  2  '' 

^Tj,k  Srj_i,k  ~  hr  ^  h  i^Jkrj—\  “i"  ^ 

(?fc)T‘('^*»-j  +  HkTj-Sr„k)',  (169) 
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if  j  is  even, 


Srj,k  Srj_^,k  ^*"2  ^  ^rj,k'Hlcj-j')  X 

For  each  rj,  it  requires  0  iNg{rj))  operations  to  solve  each  of  the  linear  systems 
(169),  (170).  Since  there  are  steps  over  [0,n7]  in  the  ODE  solve,  and  since  n^, 
ne  are  both  proportional  to  (see  (163),  (166),  and  Remark  2.1),  the  total  cost 

for  solving  the  Riccati  equation  is  0{N^). 


4.3  Evaluating  the  Operators  El-,  Er 


In  this  subsection,  we  discretize  the  linear  equation  (156)  and  evaluate  El,  Er 
on  the  grid  [pm]-  There  are  several  ways  to  evaluate  the  matrices  El,  Er;  here, 
we  present  one  of  them.  For  convenience  discussion,  we  first  scale  the  linear 
system  (156)  using  formulae  (103),  (104):  multiplying  (156)  from  the  left  by 
El{uj,  0;  Py),  and  from  the  right  by  Er{'dd,  0;  Px),  we  obtain 

ZTT  -  _ 

Jo  +  Sr,k-Hkr)-{Sq)^  X 

{Hkr-Sr,k  +  Jkr)-ER{r,  0;  Px)-r-dr 
=  EL(ro,0;Py)-(dS')o,,fc-Ei?(ro,0;Pj,).  (171) 

The  integration  over  [0,  cu]  in  (156)  is  discretized  using  the  trapezoidal  rule  on 
the  equispaced  grids  {pm}- 


•  Top 

ZTT  ^ 

—A:  hp  y  )  pjpLipj,  0,  P-u)'iJkp^  T  S p-^k'Hkpj')  x 

i=o 

i^^)pj'iEkpj'Spj^k  +  Jkpj)‘ER(^Pj,0',  P r) 

=  EL{r^,0]Py)-{SS)mxER{w,0:Px). 


where  pj  is  defined  as 


Pj  = 


Pj  J  T  2,  .  .  .  ,  Hp  1, 

Pj/2  i  =  0,np. 


(172) 


(173) 


Although  the  two  matrices  EL{r,0',  Py),  ER{r,0',  Px)  are  required  at  the  coarser 
grid  r  =  pm,  they  will  be  first  obtained  on  the  finer  grid  r  =  rj  in  order  to  main¬ 
tain  an  accuracy  comparable  to  that  in  which  the  scattering  matrix  S  is  obtained. 
We  use  formulae  (110),  (111)  to  recursively  compute  EL{r,  0;  Py),  ER{r,  0;  Px)  at 
=  '’'j,  i  =  T  2,  —  The  recursion  starts  with  the  initial  values 


f;l(0,0;P,)  =  /, 

Ej,(0,0;P,)  =  /. 
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(174) 

(175) 


For  j  =  1,2,..., Ur,  the  matrix  £i,(rj, 0;  Py)  is  updated  from  0;  Py)  via 

the  formula 

the  matrix  ER{rj,  0;  P^)  is  updated  from  0;  Px)  via  the  formula 

£;„(t-„0;P.)=  (/+hp,(r,)^  (/- jP.(ri-x))  £«(>•, -.,0;PJ.  (177) 


4.4  Forming  the  Linear  Equations 

After  EL{rj,0-,  Py),  ER{rj,0-,  Px)  are  evaluated  (see  Section  4.3),  and  Srj,k  ob¬ 
tained  from  the  forward  solve  (see  Section  4.2),  we  write  explicitly  the  linear 
equations  (172)  in  terms  of  these  matrices,  so  that  we  can  see  more  clearly  the 
structures  of  (172). 

For  j  =  1, 2, . . .  ,np,  we  denote  by  (see  Section  4.1.1) 

Nj  =  neipj)  (178) 


the  number  of  equispaced  points  on  the  circle  r  =  pj,  which  is  also  the  dimension 
of  the  matrices  Sp^^ki  EL{pj,0;  Py),  and  ER{pj,0‘,  Px)-  Given  j  and  for  0  <  m  < 
Nj,  we  further  denote  by 

_  2m; 

Nj 

the  cizimuthal  values  of  the  equispaced  points  on  the  circle  r  =  pj.  Fineilly,  for  a 
smooth  function  g,  we  denote  by  Qj^m.  the  values  of  g  at  the  equispaced  points  on 


the  circle.  We  will  require  the  matrix  of  dimension  Nn^ 

dS  =  (55)^, fc,  (180) 

and  the  three  matrices  of  dimension  Nj 

B,  =  EL{pj,^-,Py)iJkp,+Sp,,kHkp,)iN}F),  (181) 

Aj  =  {N}F-^)iJkp,AHkp,Sp„k)-ER{pj,^-,Px).  (182) 

{8q)j  =  diag{  {6q)jfl,  {Sq)j^i, ...,  {Sq)j,Nj  },  (183) 

i  i 

where  N^  F  and  Nf  F~^  are  scaled  discrete  forward  and  backward  Fourier  trans¬ 
forms,  so  that 

(iv/p)„,„  =  (184) 

(JV/p-*)„,„  =  exp  j  ■  (185) 
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Then,  with  the  help  of  (181),  (182),  (183),  (159),  we  rewrite  the  linear  equation 
(172)  in  the  form 

w-  Py)-  ^.By{6q)yA^  •Efi(0,  P,)  =  dS,  (186) 

(with  the  summation  of  matrices  of  different  dimensions  discussed  in  Remark 
4.1).  Denoting  by  V  the  linear  space  of  real- valued  vectors  of  dimension 

Tip 

=  (187) 

3-1 

and  by  P  the  linear  space  of  complex-valued  vectors  of  dimension 

(188) 

we  observe  that  dS  €  7?.,  that  (186)  is  a  system  of  linear  equations  for  the  vector 
6q  ^T>  defined  by 


^9-  {  (('^?)l, 05  (^9)1,1,  (d9)i,iVi-l), 

((^?)2,o,  (^9)2,15  ■  •  •  5  (<^9)2,Ar2-l)5 

.  .  .  , 

((^9)np,0,(<59)np,l,.-.,(d9)„p,;v„p-l)  }^,  (189) 

and  that  (186)  defines  a  linear  operator  A  :  V  P,  such  that  (172)  can  be 
reformulated  as 

A  (Sq)  =  dS.  (190) 

Remark  4.2  It  is  easy  to  see  that  the  procedures  for  obtaining  the  matrices 
{  J  ~  1,2,...  ,np  },  as  well  as  the  application  of  A  to  a  vector^  cost 

0{N^)  arithmetic  operations. 

Remark  4.3  We  define  the  inner  product  in  the  range  space  IZ  by  the  formula 

Ns 

{u,v)  =  J2uj-Vj  (191) 

i=i 

so  that  the  induced  norm  for  P  is  the  standard  norm.  Since  vectors  in  V  of 
the  form  (189)  represent  functions  of  L^{D{zu))  in  polar  coordinates,  and  since 
the  inner  product  of  a  pair  of  functions  in  L^{D{w))  defined  by  the  formula 

{f,9)  =  J^  /(r,6')-5f(r,6')d6') -r-dr,  (192) 
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induces  the  norm,  we  accordingly  define  the  discrete  version  of  (192)  by  the 
formula 

Up  /  Y 

(x,  y)r  =  hp  ^2  Pj  ( 

j=0  m=0 

as  the  inner  product  for  the  vectors  x,  y  ^V.  We  refer  to  the  norm  of  V  induced 
by  (192)  as  the  U  norm. 


(193) 


4.5  Least-square  Solution  of  the  Linearized  Equation 

In  order  to  solve  the  least-squares  problem  (190),  the  conjugate-gradient  method 
was  employed  (to  the  normal  equation  of  (190))  because  it  is  the  least  expansive 
among  several  standard  methods,  including  QR  decomposation,  Gram-Schmidt 
orthogonalization.  The  application  of  the  conjugate-gradient  method,  in  this 
case,  is  tedious  and  straightforward  with  one  exception.  The  inner  product  (193) 
must  be  used  (see  Remark  4.3)  to  obtain  the  least-squares  solution  in  V,  namely, 
the  minimization  of  the  solution  in  the  norm  of  X^,  as  well  as  the  minimization 
of  the  residual  in  the  norm  of  'll.  More  specifically,  this  means  that  the  inner 
product  in  a  standard  conjugate-gradient  method  must  be  replaced  with  (193). 
It  also  means  that  the  adjoint  operator  of  A.  (see  (190))  must  be  redefined  as 
follows.  We  denote  by  4  a  matrix  of  dimension  Ns  x  Nq  whose  i-th  row,  is 
such  that 

([oW]'.  y)^  =  M(s,a  (194) 

for  all  y  G  X>.  Obviously,  the  Hermitian  of  A,  namely,  the  linear  operator  A*  : 
7?.  hh.  P,  is  the  adjoint  operator  of  A,  with  respect  to  the  inner  product  (191)  in 
TZ,  and  to  the  inner  product  (193)  in  P.  Let  8q  =  A*u,  where  u  €  TZin  the  form 
of  an  Ns  X  Ns  matrix,  and  6q  in  the  form  (189).  It  is  a  tedious  but  straightforward 
manipulation  to  show  that 


ik^ 


Ni 


{Sq)j,m  =  ^  E  m;  P,)-u*-£i(0,  P,)  (195) 

n,l=l 

where  {T]m.,n  denotes  the  (m,  n)-th  entry  of  the  matrix  T. 

Remark  4.4  It  is  easy  to  see  from  (195)  that  the  application  of  A*  to  a  vector 
cost  0{N2)  arithmetic  operations  (see  also  Remark  f.2). 

Remark  4.5  Since  only  an  approximate  solution  of  the  least-squares  is  required 
(see  Remark  3.  If),  the  conjugate-gradient  iteration  is  usually  terminated  at  n-th 
step  in  our  numerical  experiments  when  the  ratio  of  norms  of  the  last  and  the 
initial  residuals 

a  =  M  (196) 


llroll 


is  about  10 
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5  Numerical  Results  and  Discussions 


FORTRAN  programs  have  been  written  implementing  the  procedures  described 
in  the  preceding  section.  In  this  section,  to  illustrate  the  performance  of  the  algo¬ 
rithm,  we  present  several  numerical  examples  for  the  inversion  of  the  Helmholtz 
equation  in  two  dimensions.  Remarks  will  be  made,  at  the  beginning  and  the  end 
of  this  section,  to  discuss  some  technical  details  of  the  numerical  experiments. 


5.1  The  Up- recursion  and  the  Complexity 

In  our  numerical  implementations,  only  an  approximate  is  sought  to  start 
the  recursive  procedure  described  in  Section  3.6).  We  find  consistently  in  our 
numerical  experiments  that  the  lowest  frequency  ki  can  be  chosen  such  that  its 
corresponding  wavelength  is  the  about  size  of  the  scatterer. 

Remark  5.1  Our  numerical  experiments  further  show  that  frequently  kj  (see 
Figure  7)  can  he  chosen  such  that  the  size  of  the  scatterer  is  about  j  wavelengths. 
For  instance,  we  may  set 

kj  =  j,  (197) 

for  a  scatterer,  whose  function  q  is  not  large,  inside  a  disk  of  diameter  27r. 

Assuming  a  finite  number  of  iterations  are  required  in  the  conjugate  gradient 
method,  for  the  linear  system  (157)  needs  not  be  solved  accurately  according  to 
Remark  3.14,  we  observe  that  the  least-squares  solution  of  (157)  can  be  approxi¬ 
mately  obtained  at  a  cost  of  0{N^)  arithmetic  operations  (see  Remarks  4.2,  4.4). 
Then,  the  inversion  algorithm  requires  0{N^)  operations  since  there  are  about 
frequencies  employed  in  the  recursion  (see  Remark  5.1). 

5.2  The  Forward  Modeling 

The  scattering  data  (see  Section  2.2)  are  obtained  by  numerical  solution  of  the 
forward  scattering  problem — the  initial  value  problem  of  the  Riccati  equations 
(see  (34)  and  (26)).  In  our  numerical  computation,  we  assume  the  scatterer  q  is 
nonzero  in  a  disk  of  radius  zo  =  tt. 

We  used  both  the  standard  fourth  order  Runge-Kutta  method  and  the  second 
order  implicit  scheme  described  in  Section  4.2  for  the  numerical  solution  of  the 
ordinary  differential  equation  (34).  The  numerically  obtained  solutions  were 
compared  with  the  exact  ones  when  they  were  analytically  available.  We  found 
that  both  methods  converged  at  rates  as  expected,  with  the  former  being  more 
efficient  if  the  accuracy  required  was  higher  that  10“^.  For  general  scatterers 
to  which  exact  solutions  of  the  scattering  matrix  were  not  available,  the  rates 
of  convergence  of  the  two  methods  were  verified  numerically.  In  the  numerical 
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reconstructions  presented  below,  the  scattering  data  used  were  all  obtained  with 
an  accuracy  10"^  to  10“'^. 

5.3  Numerical  Examples 

A  large  number  of  numerical  experiments  have  been  made  in  which  several  types 
of  the  scatterers  have  been  reconstructed.  The  reconstructions  of  five  types  of 
scatterers  are  presented  here.  The  biggest  problems  that  have  been  tested  are 
about  forty  wavelengths,  which  require  two  to  three  hours  CPU  time  on  a  c-90 
Cray  computer. 

Example  1:  Reconstruct  a  scatterer  defined  by 

q^{x,  y)  =  0.15-(1  -  -  0.5-  “  v") 

(198) 

inside  the  disk  T>(7r);  see  Figure  8  for  surface  and  contour  plots  of  the  scatterer 
function.  Nine  frequencies  were  used  in  the  reconstruction,  corresponding  to  wave 
numbers  k  =■  1, 2,  •  ■  • ,  9.  The  inversion  algorithm  reconstructed  it  accurately  (the 


k 

1 

2 

3 

4 

5 

6 

7 

8 

9 

62 

0.57 

0.41 

0.16 

3.1-2 

6.4E-3 

2.2E-3 

1.2E-3 

7.9E-4 

5.6E-4 

Table  1:  Error  of  Reconstruction  at  9  Frequencies,  Example  1. 

reconstructed  function  will  not  be  plotted  against  the  exact  since  the  error  is  so 
small  that  it  is  invisible  in  the  plot).  The  procedure  cost  122  seconds  CPU  time 
on  a  Cray  C-90  computer;  see  Table  1  for  the  error  of  the  reconstruction. 

Example  1.1:  To  test  the  stability  of  the  algorithm,  we  reconstruct  in  this 
example  the  scatterer  qi  but  with  noisy  data.  Noise  is  added  to  the  scattering 
data  used  in  Example  1  by  truncating  each  number  in  the  data  to  certain  digits. 
For  instance,  truncating  the  number  0.129876  to  two  digits  yields  0.12,  and  the 
perturbation  (or  noise)  incurred  here  by  the  truncation  is  about  1%.  Three  tests 
were  made  here  corresponding  to  truncations  of  the  scattering  data  to  = 
3,  2,  1  digits.  The  resulting  errors  in  the  inversion  are  listed  in  Table  2. 


Nd 

3 

2 

1 

62 

1.4E-3 

Table  2:  Error  of  Reconstruction  with  Noisy  Data,  Example  1.1. 
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Figure  9  shows  the  surface  and  contour  plots  of  the  reconstruction  with  the 
scattering  data  truncated  to  one  digit,  namely,  the  noisy  data  used  here  only 
have  one-digit  accuracy.  This  is  a  quite  severe  test  to  an  algorithm.  This  time, 
not  only  our  computation  didn’t  blow  up,  it  actually  reconstructed  the  scatterer 
with  a  11%  error.  Stability  tests  were  also  performed  to  other  scatterers  whose 
reconstructions  are  presented  in  this  paper,  the  results  being  similar. 

Example  2:  Reconstruct  a  scatterer  defined  by 

92(3;,  y)  =  0.2  {1  cos(ll-a;)  +  sin(ll-7/)}  ,  (199) 

inside  Z)(7r);  see  Figure  10  for  surface  and  contour  plots  of  the  scatterer.  This  is  a 
quite  oscillatory  function.  Inside  the  disk,  there  are  about  160  peaks  and  valleys 
representing  a  highly  rugged  index  of  refraction.  The  computation  simulates  an 
acoustic  experiment  in  which  the  background  speed  of  sound  is  that  of  water; 
the  scatterer  is  20.87  centimeters  in  diameter  (about  the  size  of  a  human  head). 
Nine  frequencies,  /  =  7, 14, 21,  •  ■  • ,  63  kHz,  were  used  in  the  reconstruction, 
corresponding  to  wave  numbers  ^  =  1, 2,  ■  •  • ,  9.  At  /  =  63  kHz,  108  transducers 
were  required  around  the  scatterer.  The  procedure  cost  122  seconds  CPU  time 
on  a  Cray  C-90  computer;  see  Table  3  for  the  error  of  the  reconstruction. 


k 

1 

2 

3 

4 

5 

6 

7 

8 

9 

^2 

0.53 

0.56 

0.56 

0.55 

0.25 

8.2E-2 

3.6E-2 

1.9E-2 

1.2E-2 

Table  3;  Error  of  Reconstruction  at  9  Frequencies,  Example  2. 

Because  of  the  complicated  structure  of  the  scatterer,  the  reconstructed  scatterer 
is  plotted  against  the  exact  scatterer  first  horizontally  across  the  diameter  of  the 
disk  T)(7r),  then  across  concentric  circles  of  various  radii  between  0  and  tt  of  the 
disk;  see  Figures  11,  12. 


Example  3:  Reconstruct  a  scatterer  defined  in  i9(7r)  by 


<l3{x,y) 


'  qi{x  10.8, y  10.8) 
<  -0.5 
^  0 


if  r  <  2.6, 
if  r  G  [2.6, 2.9), 
if  r  >  2.9; 


(200) 


see  Figure  13  for  surface  and  contour  plots  of  the  function.  This  scatterer  is 
difficult  to  reconstruct  for  two  reasons. 


(1)  Across  the  two  circles  r  =  2.6  and  r  =  2.9,  the  function  is  discontinuous. 
The  value  of  the  function  changes  sharply  to  —0.5  in  the  narrow  annulus. 
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(2)  If  the  background  speed  of  sound  is  that  of  water,  then  the  material  in 
the  narrow  band  2.6  <  r  <  2.9  has  speed  of  sound  1.4  times  as  large  as  that 
of  water.  As  a  result,  this  high-speed  region  with  sharp  boundaries  blocks  the 
passage  of  the  probing  sound  waves  to  the  inside  of  the  structure,  making  it  hard 
to  reconstruct  the  smooth  part  of  the  scatterer  in  the  middle  of  the  object. 

This  example  could  be  regarded  as  a  model  problem  for  ultrasound  tomography 
of  a  human  head,  where  the  skull  is  represented  by  the  thin  layer  of  denser 
material  in  the  region  2.6  <  r  <  2.9.  If  the  actual  object  is  20.87  centimeters  in 
diameter,  the  frequencies  used  were  /  =  7, 14, 21,  •  •  • ,  84  kHz,  corresponding  to 
wave  numbers  A;  =  1, 2,  •  •  • ,  12.  At  /  =  84  kHz,  128  transducers  were  used  around 
the  scatterer.  The  CPU  time  required  for  the  procedure  was  263  seconds  on  a 
Cray  C-90  computer.  The  errors  of  the  reconstruction  at  the  12  frequencies 


k 

1 

2 

3 

4 

5 

6 

62 

0.576 

0.510 

0.367 

0.260 

0.231 

0.220 

k 

7 

8 

9 

10 

11 

12 

62 

0.197 

0.164 

0.146 

0.141 

0.138 

0.136 

Table  4:  Error  of  Reconstruction  at  12  Frequencies,  Example  3. 

are  listed  in  Table  3.  Figure  14  shows  the  surface  and  contour  plots  of  the 
reconstructed  scatterer,  whereas  Figure  15  shows  the  reconstruction  horizontally 
across  the  diameter  of  the  scatterer. 

An  examination  of  the  plots  show  that  the  error  of  the  reconstructions  occurs 
largely  around  the  discontinuities,  while  the  smooth  part  is  recovered  more  ac¬ 
curately. 

Example  4:  Reconstruct  a  cylindriceilly  symmetric  scatterer  defined  by 

^^(r)  =  0.3  •  [  0.55cos(2r)  —  0.44sin(4r)  -|-  0.23sin(6r)  -f-  0.3cos(8r)  ],  (201) 

for  0  <  r  <  TT.  Nine  frequencies  were  used  in  the  reconstruction,  corresponding 
to  wave  numbers  k  =  1,2, 3, ••■,9.  The  purpose  of  this  example  is  to  show 
how  the  Uncertainty  Principle  works  by  illustrating  the  process  of  convergence 
and  the  distributions  of  error  in  reconstructions  at  these  nine  frequencies.  An 
examination  of  the  plots  of  the  reconstructions  and  error  functions  shows  that 

(1)  The  largest  error  of  the  reconstructions  occurs  near  the  two  points  r  =  0, 
r  =  TT  where  the  scatterer  is  not  smooth. 

(2)  The  error  function  q  —  qk  is  more  or  less  evenly  distributed  away  from  the 
end  points  r  =  0,  r  =  tt. 
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(3)  The  higher  the  wave  number  k  is,  the  more  oscillatory  the  error  function 
becomes. 


Example  5:  Reconstruct  a  cylindrically  symmetric  scatterer  in  the  disk  D{'k) 
defined  as  follows.  Give  tq  >  0,  the  scatterer  function  is  given  by  the  formula 


q5{r)  =  0.4  • 


1  —  cos 


(202) 


for  To  <  r  <  -K.  Inside  the  disk  /^(t’o),  the  scatterer  function  q  is  not  defined; 
on  the  circle  r  =  ro,  the  Dirichlet  boundary  condition  is  imposed  on  the  total 
field  <f).  Therefore,  the  circle  |x|  =  ro  is  a  totally  reflecting  surface,  and  the 
medium  inside  D{ro)  is  “sound-soft”  (see,  for  example,  [1],  page  319-321  for  the 
treatment  of  exterior  Dirichlet  problem).  We  choose  ro  =  0.8  in  the  numerical 
computation.  Fifteen  frequencies  were  used  in  the  reconstruction,  corresponding 
to  wave  numbers  ^  =  1, 2, 3,  •  ■  • ,  12. 

The  inversion  algorithm  doesn’t  know  the  fact  that  part  of  the  scattering  is 
due  to  a  reflecting  surface.  Rather,  it  assumes  that  the  scattering  occurs  as  a 
result  of  a  distribution  of  an  unknow  index  of  refraction  everywhere  continuous: 
on  the  disk  D(7r),  across  the  circle  r  =  ro,  and  inside  the  disk  Z)(ro).  The 
purpose  of  this  example  is  to  see  whether  the  inversion  algorithm  will  be  stable 
and  convergent,  and,  if  so,  what  this  reconstructed  continuous  function  q  will  be, 
particularly  across  the  reflecting  surface  and  inside  the  disk  D{ro). 

Our  numerical  experiments  show  that  the  reconstruction  is  indeed  stable  and 
convergent,  see  Figure  18  for  the  plots  of  the  reconstructed  q^  depicted  against 
the  exact  q  at  the  fifteen  frequencies.  In  the  pictures,  the  exact  q  is  set  zero  in 
the  interval  [0,  ro]  to  indicate  that  it  is  not  defined  inside  Zl(ro).  The  errors 
given  in  Figure  18  are  measured  in  the  interval  [ro,7r].  Two  comments  can  be 
made  about  the  reconstruction: 


(1)  Inside  the  totally  reflecting  disk,  there  are  regions  where  the  values  of  the 
reconstructed  qk  are  less  than  —1.  After  the  reconstruction  starts  to  converge, 
namely,  when  k  is  greater  than  4,  ^^(r)  decreases  as  r  moves  across  the  reflecting 
disk,  and  becomes  less  than  —1  the  moment  r  moves  into  the  reflecting  circle 
r  =  tq.  In  other  words,  the  inversion  algorithm  sees  the  reflecting  disk  as  a 
region  where  the  speed  of  sound  of  the  medium  is  imaginary  according  to  the 
formula 

k^l  -f  9(r)  =  (203) 

where  u;  >  0  is  the  angular  frequency,  c(r)  >  0  is  the  speed  of  sound.  Therefore,  to 
the  full  aperture  experiment  where  the  scattering  data  are  collected,  a  reflecting 
disk  looks  the  same  (see  Definition  3.8)  as  a  layer  of  non-propagating  medium. 

(2)  Outside  the  totally  reflecting  disk,  the  smooth  part  of  the  prescribed  scat¬ 
terer  is  essentially  recovered:  the  reflecting  surface  is  regarded  as  a  discontinuity 
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of  some  type,  which  causes  the  Gibbs  phenomenon  around  it;  the  scatterer  func¬ 
tion  (202)  is  reconstructed  more  accurately  away  from  this  singular  point. 


5.4  Discussions  and  Conclusions 

The  following  technical  details  of  the  numerical  implementation  appeax  to  be 
worthy  mentioning. 

1.  The  convergence  of  recursive  linearization  procedure  depends,  of  course,  on 
the  step  size  of  the  frequency  k.  We  find  in  our  numerical  experiments  that  when 
the  scatterer  function  is  not  large  (for  example,  when  —  0.8  <  ?  <  1),  convergence 
is  usually  guaranteed  with  the  step  size  given  by  (197).  In  general,  smaller  step 
sizes  are  required  when  the  scatterer  function  is  very  large  or  very  closed  to  —  1. 

On  the  other  hand,  larger  step  sizes  of  k  can  generally  be  used  at  higher  frequen¬ 
cies  k  without  affecting  the  convergence.  This  is  so  because  at  a  relatively  high 
frequency  k  =  a  where  the  dominant  lower-frequency  components  of  the  scatterer 
have  been  recovered,  q  —  Qa  is  small.  Therefore,  the  perturbation  6q  =  qb  —  qa 
will  be  small  for  a  relatively  large  step  size  of  8k  =  h  —  a. 

2.  The  stability  of  the  algorithm  is  not  sensitive  to  the  step  size  of  the 
frequency  k.  It  is  largely  controlled  by  the  way  the  ill-posed  linear  system  (157) 
is  solved.  Numerical  experiments  show  that  the  up-recursion  in  frequency  k  is 
usually  unstable  when  the  least-squares  solution  of  (157)  is  obtained  in  such  a 
precision  that  a  (see  (196))  is  smaller  than  10“®.  With  scattering  data  accurate 
to  three  digits,  we  find  that  a  =  10~^  (see  Remark  4.5)  is  suitable  to  inversion  in 
a  varieties  of  cases. 

The  following  discussions  are  about  the  models  of  forward  scattering. 

3.  The  use  of  cylindrical  geometry  to  introduce  the  scattering  matrix,  and 
subsequently,  to  obtain  the  Riccati  equation,  is  a  convenient  but  not  the  only 
approach  to  the  forward  modeling.  Scattering  matrix  associated  with  straight- 
line  geometry,  and  its  Riccati  equation,  for  example,  are  introduced  in  [10]. 

4.  Forward  models  other  than  the  Riccati  equations  can  be  used  to  recursively 
linearize  the  inverse  problem.  The  Lippmann- Schwinger  equation  seems  a  better 
candidate  for  the  forward  modeling  because  there  is  a  recursive  procedure  that 
solves  accurately  the  forward  problem  in  0{N^)  operations  at  a  frequency.  Thus, 
the  multi-frequency  inversion  costs  about  0{N^)  operations.  The  linearization 
procedure  based  on  the  Lippmann-Schwinger  equation  is  described  in  [12].  Its 
implementation  and  numerical  results  will  be  reported  on  a  later  date. 

The  following  discussions  are  about  extentions  of  the  algorithm  to  three  dimen¬ 
sions  and  to  other  types  of  scattering  problems. 
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5.  The  direct  extention  of  the  recursive  linearization  procedure  to  three  di¬ 
mensions  is  straightforward,  but  is  difficult  to  be  implemented  numerically,  for 
the  procedure  requires  0{N^)  operations.  The  high  computational  cost  in  three 
dimensions  is  primarily  a  result  of  the  so-called  data  redundancy:  at  a  frequency 
A:  >  0,  the  full-aperture  scattering  data  depend  on  four  independent  parameters 
whereas  the  scatterer  q  to  be  recovered  is  a  function  of  three  spatial  parameters. 
There  are  several  ways  to  reduce  the  cost  to  0{N^)',  each  of  which  uses  part  of 
the  full-aperture  scattering  data. 

6.  The  recursive  procedure  can  be  applied  to  measurements  of  limited  aper¬ 
ture  (see  [12]  for  more  details).  In  two  dimensions,  this  only  allows  us  to  recover 
partial  information  of  the  scatterer  within  the  Fourier  aperture  D{2k)'^ .  In  three 
dimensions,  measurements  of  limited  aperture  may  provide  the  same  amount  of 
information  about  the  scatterer  as  do  the  full-aperture  measurements. 

7.  The  scheme  can  be  used  to  solve  inverse  scattering  problems  of  more 
complicated  equations  describing  more  realistic  processes  of  acoustic,  elastic,  or 
electromagnetic  scattering  in  which  the  Heisenberg’s  Uncertainty  Principle  holds. 
The  exact  formulation  of  the  uncertainty  principle  may  differ  in  specific  environ¬ 
ments,  but  it  is  certain  that  everywhere  in  the  realm  of  wave  phenomena  an 
incident  wave  of  a  lower  frequency  interacts  weakly  with  the  Fourier  modes  of 
the  scatterer  of  higher  frequencies,  and  such  an  interaction  produces  a  weaker 
scattered  field. 

8.  Heisenberg’s  Uncertainty  Principle  is  also  valid  in  the  case  of  obstacle  scat¬ 
tering.  There,  lower-frequency  incident  fields  interacts  weakly  with  the  higher- 
frequency  roughness  of  the  surface  of  the  scattering  obstacle.  Therefore,  the 
inverse  obstacle  scattering  problem  can  be  recursively  linearized  by  the  same 
mechanism  introduced  here  in  this  paper. 
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Figure  9:  Surface  and  Contour  Views  of  Inversion  with  iVj  =  1,  Example  1.1. 
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Figure  10:  Surface  and  Contour  Views  of  Scatterer  52,  Example  2. 
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Figure  11:  Reconstructed  vs  Exact  52  on  Diameter,  k  =  1, 2, ...  9,  Example  2. 
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Figure  13;  Surface  and  Contour  Views  of  Scatterer  qz.  Example  3. 
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Across  Diameter.  K=10.0  err=0.Ul  Across  Dlomeler,  K=11.0  err=0.I38  Across  Dlonieter,  K=12,0  err^O.lSS 


Figure  15:  Reconstructed  vs  Exact  on  Diameter  at  12  Frequencies,  Example  3. 


49 


